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Abstract 

We  consider  an  Eulerian-Lagrangian  localized  adjoint  method  (ELLAM)  applied 
to  nonlinear  model  equations  governing  solute  transport  and  sorption  in  porous 
media.  Solute  transport  in  the  aqueous  phase  is  modeled  by  standard  advection  and 
hydrodynamic  dispersion  processes,  while  sorption  is  modeled  with  a  nonlinear  local 
equilibrium  model.  We  present  our  implementation  of  finite  volume  ELLAM  (FV- 
ELLAM)  and  finite  element  (FE-ELLAM)  discretizations  to  the  reactive  transport 
model  and  evaluate  their  performance  for  several  test  problems  containing  self- 
sharpening  fronts. 
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Notation 

Roman  Letters 

C  aqueous  phase  solute  concentration 

C  trial  solution 

Cb  inflow  boundary  value  for  aqueous  phase  solute  concentration 

Ce  aqueous  phase  equilibrium  solute  concentration 

Cn+1’m  nonlinear  solver’s  current  guess  for  concentration  at  tn+1 

C*  boundary  value  at  —  oo  for  traveling  wave  example 

C*  boundary  value  at  oo  for  traveling  wave  example 

Cl  Riemann  problem  left  state 

CR  Riemann  problem  right  state 

D  hydrodynamic  dispersion  coefficient 

K f  Freundlich  sorption  capacity  coefficient 

M  normalized  total  concentration,  C  +  tp(C) 

A4h  discrete  spatial  mesh 

Rf  retardation  factor,  1  + 
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T  extent  of  temporal  domain 

W  FV-ELLAM  approximate  test  function 

X  parameterization  for  tracking  along  solution  characteristics  (see 

eqns  (27)  and  (28)) 

Xa  parameterization  for  tracking  along  adjoint  characteristics  (see 

eqns  (25)  and  (26)) 

Cfr  wave  speed  cutoff  for  FT  tracking 

d\-e  front  width  in  traveling  wave  example 

/  continuous  solute  sources  and  sinks 

n  outward  unit  normal  on  dtt 

ne  +  1  number  of  nodes  in  Xih 

rif  exponent  in  Freundlich  sorption  isotherm 

Tiq  number  of  integration  points 

nst  number  of  tracking  substeps 

qb  total  flux  boundary  value  on  inflow  boundary 

rq  numerical  integration  weight 

t  time  coordinate 

tn,k  tracking  time  level,  tn,k  =  tn  +  kAr 

v  mean  pore  velocity 

w  test  function 

x  space  coordinate 

xl  left  end  of  spatial  domain 

Xr  right  end  of  spatial  domain 

Xi+ 1/2  FV-ELLAM  cell  boundary,  ( Xi  +  Xi+ 1)/2 

x™+1  intersection  of  shock  at  time  level  n  +  1  in  Riemann  example 

xq  numerical  integration  point 


3 


Greek  Letters 


A  tn+1 


At 


A  Xi 


Ax 

Ax 

Ax 

r0 

r7 

n 


*+1/2 


min 


Vi 


a 

dtl 

5fr 

e 


V 

e 

A 

Aa 

Pb 

a 


<P 

LO  e 


Uf 


time  step,  tn+1  —  tn. 

tracking  time  step 

spatial  increment,  xi+\/2  —  Xi-1/2 

spatial  increment,  xl+\  —  xt 

support  parameter  for  Wt ,  Ax  =  Ax/NS 

approximate  minimum  front  resolution  required 

outflow  spatial  boundary 

inflow  spatial  boundary 

spatial  domain 

spatial  interval  associated  with  node  % 

wave  speed  for  traveling  wave  example 

boundary  of  spatial  domain 

discretization  parameter  for  FT  tracking 

front  width  parameter  for  traveling  wave  example 

traveling  wave  coordinate,  r]  =  x  —  at 

porosity 

characteristic  speed,  v/Rf 

adjoint  equation  characteristic  speed,  (vC)/[C  +  <p(C)\ 
bulk  density  of  the  solid  phase 
shock  speed 

normalized  isotherm,  \pbu>e(C)\/ 6 
piecewise-linear  Lagrangian  shape  function 
solid  phase  equilibrium  solute  mass  fraction 
solid  phase  solute  mass  fraction 
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Subscripts  and  Superscripts 


a 

i 

J 

m 

n 


Abbreviations 

BE 

BE-S 

BL 

Cr 

ELLAM 

FE-ELLAM 

FT 

FV-ELLAM 

MMOC 

NS 

NT 

Pe 

RK2 

RK2-S 


quantity  associated  with  adjoint  characteristics  (subscript) 
nodal  value  (subscript) 
nodal  value  (subscript) 

nonlinear  solver  iteration  level  (superscript) 
time  level  (superscript) 
forward-tracked  quantity  (superscript) 
backward-tracked  quantity  (superscript) 


backward  Euler  time  discretization 
tracking  strategy  using  BE  with  At  =  Atn+1 
bilinear  interpolation  in  space  and  time 
Courant  number 

Eulerian-Lagrangian  localized  adjoint  method 

finite  element  ELLAM 

front-tracking  method 

finite  volume  ELLAM 

modified  method  of  characteristics 

FV-ELLAM  parameter  for  approximate  test  function,  W 
number  of  composite  trapezoidal  rule  intervals  (in  time)  along 
inflow  boundary 
mesh  Peclet  number 

second-order  explicit  Runge-Kutta  time  discretization 
tracking  strategy  using  RK2  with  At  =  Atn+1 
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RK2-BL 


tracking  strategy  combining  RK2  with  BL 
RK2-FT  tracking  strategy  combining  RK2  with  FT 
SSIP  strategic  spatial  integration  point 

STIP  strategic  temporal  integration  point 

1  Introduction 

Much  effort  has  been  devoted  to  the  numerical  simulation  of  contaminant 
transport  processes  in  the  subsurface  over  the  last  few  decades  [5,  28].  Despite 
significant  advancements,  accurate  and  efficient  simulation  remains  a  challenge 
in  many  cases,  particularly  for  advective-dominated  problems  involving  nonlin¬ 
ear  chemical  reactions  and  mass  transfer  [16,  20].  Characteristic-based  meth¬ 
ods  such  as  Lagrange-Galerkin  discretizations  [4,  29],  the  modified  method  of 
characteristics  (MMOC)  [13],  and  ELL  AMs  [9]  have  been  applied  to  a  wide 
range  of  transport  problems.  These  methods  typically  combine  a  Lagrangian 
approach  for  advection  with  Eulerian  discretizations  for  other  transport  pro¬ 
cesses  such  as  physical  dispersion.  Since  they  rely  on  a  Lagrangian  framework 
for  advection,  characteristic-based  methods  are  often  able  to  provide  sharp  res¬ 
olution  of  fronts  on  relatively  coarse  grids  while  avoiding  stability  restrictions 
on  the  Courant  number  commonly  found  in  Eulerian  methods  [14],  Among 
characteristic-based  methods,  ELLAM  approaches  have  the  additional  advan¬ 
tages  that  they  provide  mass  conservation  and  incorporate  boundary  condi¬ 
tions  in  a  systematic  way  [31]. 

A  general  review  of  characteristic-based  methods  and  ELLAM  approaches 
in  particular  is  beyond  the  scope  of  this  work.  We  refer  the  reader  instead 
to  Russell  and  Celia  [31]  and  Ewing  and  Wang  [14],  In  brief,  ELLAM  dis- 
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cretizations  for  linear  transport  are  mature.  A  number  of  practical  two-  and 
three-dimensional  problems  have  been  solved  successfully  using  ELLAM  ap¬ 
proximations  [7,  8,  18,  40],  and  an  FV-ELLAM  option  has  been  incorporated 
in  the  USGS  MOC3D  code  [32],  ELLAM  approaches  have  also  been  consid¬ 
ered  for  several  nonlinear  problems  including  two-phase  flow  via  the  Buckley- 
Leverett  equation  with  capillary  pressure  [11],  and  advective-dispersive  trans¬ 
port  with  biodegradation  [6,  35,  41].  Initial  work  has  been  reported  for  general 
compositional  formulations  of  multiphase  flow  [42]  as  well.  However,  ELLAM 
approaches  for  nonlinear  problems  are  less  advanced  than  their  linear  counter¬ 
parts.  Nonlinear  transport  problems  introduce  several  additional  complexities. 
Solution  fronts  can,  for  example,  be  self-sharpening  (shocks)  rather  than  con¬ 
tact  discontinuities,  and  reaction  terms  can  be  sensitive  to  overshoot  and  un¬ 
dershoot  in  the  solution  [27].  To  realize  the  same  level  of  success  that  has  been 
achieved  for  linear  problems  with  nonlinear  reactive  transport,  more  work  is 
needed  to  develop  approaches  within  the  ELLAM  framework  that  can  effec¬ 
tively  maintain  accurate  resolution  for  large  time  steps  and  mass  conservation 
in  the  presence  of  these  additional  challenges  [8]. 

Previous  ELLAM  approaches  for  nonlinear  reactive  transport  have  consid¬ 
ered  contaminant  biodegradation  modeled  by  Monod  kinetics  in  both  fully 
coupled  [35,  41]  and  operator-split  frameworks  [6].  Within  an  operator  split¬ 
ting  context,  ELLAM  methodologies  can  be  carried  over  directly  to  the  linear 
transport  equation.  The  fully  coupled  formulations  account  for  the  reaction 
terms  in  the  ELLAM  test  functions.  The  methods  of  Wang  et  al.  [41]  and 
Vag  et  al.  [35]  are  based  on  linearizing  the  Monod  reaction  terms  and  defining 
ELLAM  test  functions  to  increase  or  decrease  along  characteristics,  as  with 
solution  to  problems  involving  transport  and  linear  decay.  Various  techniques 
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are  then  necessary  to  account  for  the  error  arising  from  the  linearization.  How¬ 
ever,  the  basic  methods  for  tracking  information  along  characteristics  in  the 
corresponding  linear  transport  problem  can  still  be  used  [35,  41].  For  many 
nonlinear  problems  of  interest,  this  approach  will  not  be  effective. 

Several  Eulerian-Lagrangian  discretizations  including  characteristic  Galerkin 
[12]  and  Lagrange-Galerkin  methods  [4,  23,  24]  have  been  applied  to  one¬ 
dimensional  advective-dispersive  transport  with  both  equilibrium  and  non¬ 
equilibrium  sorption  in  non-divergence  form.  ELLAM  formulations  that  are 
applied  directly  to  the  divergence  form  of  the  underlying  transport  problem 
have  advantages  in  handling  of  boundary  conditions  and  ensuring  conservation 
of  mass  [31].  In  addition,  previous  Eulerian-Lagrangian  methods  accounted 
for  Lagrangian  aspects  of  their  discretization  using  characteristics  from  the 
hyperbolic  portion  of  the  original  transport  equation  [4,  23,  24],  While  tracking 
these  characteristic  curves  in  space  and  time  details  the  propagation  of  solution 
values,  characteristics  become  nonunique  in  the  presence  of  a  shock  [27].  A 
potential,  under-explored  solution  to  this  problem  is  to  track  characteristics 
for  the  ELLAM  test  functions,  which  satisfy  a  linear  adjoint  equation,  and 
so  avoid  the  difficulties  associated  with  tracking  characteristics  of  the  original 
transport  equation  when  problems  exhibit  self-sharpening  fronts. 

The  overall  goal  of  this  work  is  to  advance  ELLAM  approaches  for  the  solu¬ 
tion  of  a  common  class  of  nonlinear  transport  problems  using  approaches  that 
conserve  mass,  are  able  to  resolve  self-sharpening  fronts,  and  accommodate 
boundary  conditions  naturally.  The  specific  objectives  of  this  work  are:  (1) 
to  summarize  a  common  nonlinear  transport  model  that  poses  challenges  for 
traditional  ELLAM  approaches;  (2)  to  detail  mass  conservative  FE-ELLAM 
and  FV-ELLAM  approximation  of  the  nonlinear  problem;  (3)  to  formulate 


alternative  approaches  for  approximating  the  nonlinear  tracking  along  char¬ 
acteristics;  and  (4)  to  compare  the  various  ELLAM  approaches  for  a  range  of 
test  problems  and  spatial  and  temporal  step  sizes. 


2  Model  Formulation 


We  sought  a  nonlinear  model  problem  that  embodied  the  challenges  of  sharp- 
front  resolution,  mass  conservation,  and  the  flexible  and  accurate  accommo¬ 
dation  of  boundary  conditions,  which  were  highlighted  above  as  open  issues. 
A  common  nonlinear  model  that  meets  these  criteria  is  advective-dispersive 
transport  in  the  presence  of  nonlinear,  local-equilibrium  sorption  to  a  fixed 
solid  phase,  which  is  described  by  the  Freundlich  equilibrium  model  [15,  20]. 
This  model  is  relevant  because  it  may  be  used  as  at  least  a  first-cut  approxima¬ 
tion  for  the  transport  of  a  large  number  of  neutral  hydrophobic  solutes  through 
porous  media  that  include  soils,  sediments,  and  aquifer  materials  [2,  5].  A 
weakness  in  this  model  is  that  it  assumes  that  the  solute  achieves  equilibrium 
locally,  or  rapidly  in  comparison  to  the  rate  of  transport  through  the  system — 
an  assumption  that  may  not  hold  for  many  frequently  encountered  situations 
[5], 

Because  our  focus  is  on  methods  development,  we  examine  a  one-dimensional 
form  of  this  model  given  by 


dC  pbduf  d  (vC) 
dt  9  dt  dx 


i(»E) 

UJf 


f(x,  t )  in  x  [0,  T] 

COe(C) 


(1) 

(2) 


where  C  is  the  aqueous-phase  solute  concentration,  t  is  time,  pb  is  the  bulk 
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density  of  the  solid  phase,  6  is  the  porosity,  ujf  is  the  solid-phase  solute  mass 
fraction,  v  is  the  mean  pore  velocity,  x  is  the  spatial  coordinate,  D  is  the 
hydrodynamic  dispersion  coefficient,  /  represents  a  solute  source  or  sink,  ff  is 
the  spatial  domain,  [0,  T]  is  the  temporal  domain,  and  cce  is  the  solid-phase 
solute  mass  fraction  in  equilibrium  with  the  fluid-phase  equilibrium  solute 
concentration  Ce.  The  solid-phase  equilibrium  solute  mass  fraction  is  described 
using  the  Freundlich  equilibrium  model 

ue  =  KfCnf  (3) 

where  Kf  is  a  sorption  capacity  coefficient,  and  the  exponent  rif  is  a  measure 
of  the  sorption  intensity. 

The  initial  aqueous-phase  solute  concentration  and  solid-phase  solute  mass 
fraction  are  denoted  by  C(x,t  =  0)  =  C°(x)  and  c Of(x,t  =  0)  =  ue(C°). 
ELLAM  formulations  naturally  incorporate  a  range  of  boundary  conditions 
[31].  However,  we  restrict  ourselves  to  a  total  flux  condition  on  the  inflow 
boundary  and  zero  dispersive  flux  along  the  outflow  boundary  for  simplicity 

i^vC  —  '  n  =  Qb{x>  t)  for  x  G  T  j, 

dC 

— —  •  n  —  0  for  x  e  T0  (4) 

ox 

where  T/  U  To  =  dVL,  T/  fl  To  =  0,  and  <9h2  is  the  boundary  of  Q.  n  is  the 
outward  unit  normal  on  dQ,  with  u  ■  n  <  0  on  Tj  and  r  ■  n  >  0  on  To.  In  the 
the  one-dimensional  problems  we  consider  here,  T r  and  To  each  consist  of  a 
single  point. 

Following  traditional  ELLAM  approaches  [9,  31]  and  our  preliminary  efforts 
on  this  model  problem  [15],  we  rewrite  eqn  (1)  in  a  slightly  more  general  form 
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for  convenience: 


dM(C ) 


dt 


k  ivC 


dC\ 

Dlh)  =  ^ X ’  ^  hl  °  X  f°’  ^ 


Pb 


M  (C)  =  C+^cue(C)=C  +  ip(C) 

u 


(5) 


We  next  multiply  eqn  (5)  by  a  test  function  w(x,t),  integrate  over  fix  [0,  T] 
and  expand  to  obtain 


si  o  v  7  on 


w  vC  —  D 


dCf 

dx 


-If  (vCfr  ~  dx  dt  =  I  I  fw  dx  dt 

J  J  \  ox  ox  ox  J  J 

o  n  v  7  o  n 

After  reordering  terms,  we  have 


dx  dt 


(6) 


a\M{C)w 


n  o 


dt 


ii  /'  /'  „  dC  dw 
dtdx+  I  I  D  — — — —  dx  dt  -\- 
dx  dx 


o  n 


o  si 


o  n 


dw 


M{c)~it+vcd£ 


i 

dx  dt  +  J  J  fw  da;  dt 


9  (  r  ndC  ' 
—  vCw  —  U——w 

OX  \  OX  i 


(7) 


o  n 


Next  the  temporal  domain  [0, T]  is  divided  into  intervals  [tn,tn+1].  The  test 
function  w  is  required  to  disappear  for  t  ^  [tn,tn+1]  and  to  satisfy  the  formal 
adjoint  equation 

dw  dw 

M^m+vC  &r°  <8> 

over  [tn,tn+1].  While  other  choices  of  w  are  possible  [33],  eqn  (8)  corresponds 
to  the  standard  approach  [31].  It  is  important  to  note  that  eqn  (8)  is  nonlinear 
in  the  concentration  C  but  is  linear  in  the  test  function  w.  Assuming  that  the 
concentration  is  continuous  over  the  interval  [ tn,tn+1 }  and  applying  Green’s 
formula,  we  then  have 


da;  dt 
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(9) 


fn+ 1 

c  c  r  c)C^  rhn 

J  M[C(x,  tn+1)\w(x,  tn+1)  dx  +  J  JD-^-—dxdt 

n  tn  n 


tn+1 


tn  an 


'  r  ndC  N 
vCw  —  D——w 

,  OX  , 


tn+ 1 

+  J  J  fwdxdt 

tn  n 


n  ds  dt 


J  M[C{x,tn)}w{x,tn)dx 
n 


Eqn  (9)  consists  of  five  integrals  corresponding  to  the  total  mass  at  the  new 
time  level,  dispersion,  boundary  fluxes,  mass  at  the  old  time  level,  and  con¬ 
tributions  from  sources  and  sinks.  There  are  several  appealing  features  of  this 
weak  formulation.  Specifically,  it  is  in  a  fully  conservative  form  and  does  not 
involve  differentiating  the  isotherm  uie(C),  which  is  not  Lipschitz  continuous 
at  the  origin  for  a  Freundlich  isotherm  with  0  <  n f  <  1  [4];  the  range  of 
non-Lipschitz  continuity  occurs  routinely  in  applications  [2], 

Global  mass  conservation  for  our  ELLAM  approach  can  be  seen  by  summing 
eqn  (9)  over  all  test  functions  {wy}.  If  the  test  functions  {ny}  are  required 
to  satisfy  J2iwi(x,t)  =  1,  we  obtain  a  statement  of  mass  conservation  for  the 
domain  hi  and  the  time  interval  [tn,tn+1]: 


J  M[C(x,tn+1)]dx  +  J  J  (vC-D^-  j  -  ndsdt 

f!  t™  <90  ^  ' 

f n+1 

J  M[C(x,  tn)]  dx  +  J  J  fdxdt 


tn  n 


(10) 


Before  introducing  the  discrete  approximation  to  eqn  (9),  we  identify  the  char¬ 
acteristics  associated  with  the  original  transport  equation  and  adjoint  equa¬ 
tion,  since  these  play  an  important  role  in  the  ELLAM  approximation.  The 
hyperbolic  portion  of  eqn  (1)  has  the  characteristic  speed 
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‘-j; 

*-‘-rc 


(11) 

(12) 


where  Rf  is  the  retardation  factor  [20].  The  characteristics  are  the  solution  to 

(13) 


ClT  V 

—  =  A  [C(x,t),x,t\  =  A  (C)  = 


cLC 


The  characteristics  for  the  adjoint  eqn  (8)  are  given  instead  by 

n  nr 

—  =  A  a[C{x,t),x,t]  =  A  a[C{x,t)} 


(14) 


where  the  speed  is 


“  C  +  ip(C) 

In  general,  both  A  and  Aa  are  nonlinear  functions  of  the  concentration  C.  For 
the  Freundlich  isotherm  from  eqn  (3)  we  have 


A  = 


_|_  nfKfPb  (Jn.f-1  ’ 


and 


A„  = 


1  +  !£aiCni-1 


(16) 


When  0  <  rif  <  1.0,  which  is  the  case  we  will  focus  on  in  our  numerical 
experiments,  ^  — >  oo  as  C  \  0  and  both  A  and  Aa  — >  0.  In  addition,  |Aa|  < 
| A |  <  |u|  for  0  <  rif  <  1.0,  since  Kf  >  0  by  dehnition. 


To  provide  a  simple,  concrete  illustration  of  the  solution  and  adjoint  charac¬ 
teristic  behavior  around  a  front,  we  consider  a  Riemann  problem  with  left  and 
right  concentrations  Cl  =  0.9  and  Cr  =  0.1.  Figures  1  and  2  show  A  and  Aa  for 
v  —  1,  D  —  0.0,  and  tp(C)  =  0.50085  x  C0  '.  The  solution  is  a  right-going  shock 
moving  with  Rankine-Hugoniot  speed  a  =  v(Cr—Cl)/[M(Cr)  —  M(Cl)\.  The 
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Fig.  2.  Aa  Riemann  example 

characteristics  defined  by  eqn  (11)  are  straight  lines  in  space-time  at  which 
a  constant  concentration  value  propagates.  These  intersect  to  form  the  shock 
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seen  in  Figure  1.  On  the  other  hand,  the  adjoint  characteristics  propagate 
constant  values  of  the  test  function  w,  and  are  not  straight  since  they  depend 
on  the  solution  C.  At  the  shock,  the  adjoint  characteristic  speed  jumps  along 
with  the  solution  value.  The  characteristic  trajectory  itself  remains  continuous, 
and  the  characteristics  do  not  intersect.  In  the  presence  of  physical  dispersion, 
the  solution  characteristics  may  become  close,  but  they  will  not  cross.  Simi¬ 
larly,  the  adjoint  characteristic  speeds  will  continue  to  change  rapidly  around 
a  sharp  front  but  will  not  undergo  a  jump  in  speed  at  the  front. 


3  Solution  Approach 

Our  first  step  in  developing  a  discrete  approximation  for  the  weak  formulation, 
is  to  introduce  a  temporal  approximation  for  the  dispersion  and  source  terms  in 
eqn  (9).  Although  second-order  Runge-Kutta  methods  have  been  used  in  some 
cases  [1] ,  the  most  common  approach  is  to  use  a  backward  Euler  discretization 
in  time  so  that  the  source  and  dispersion  integrals  only  involve  values  at 
tn+i  [3xj  _  With  the  backward  Euler  approximation,  we  have  the  semi-discrete 
system 


J  M[C(x ,  tn+1)]w(x ,  tn+1 )  dx  +  A tn+1 

n 


dw 


dx 


^n+1 

+  /  /  (vCw  -  D^-wj 


n  ds  d  t 


+A tn+1  j  f(x,  tn+1)w(x,  tn+1 )  dx 
n 


J  M[C(x,tn)]w(x,tn)  dx 
n 


(17) 


where  A  tn+1  =  tn+1  —  tn. 


We  next  must  decide  on  a  representation  for  the  approximate  solution  to 
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C  and  the  test  function  w.  Our  approach  is  essentially  the  same  whether  a 
FV-ELLAM  or  FE-ELLAM  method  is  used.  In  the  following,  we  detail  both 
approaches  and  note  differences  that  exist. 

Since  the  problems  considered  here  are  one-dimensional,  the  spatial  domain 
is  simply  fl  =  [xl,xr].  We  introduce  a  discrete  mesh  M.h  for  fl  consisting  of 
nodes  or  vertices  Xi,  i  —  0, . . . ,  ne.  M.h  is  a  natural  description  of  the  compu¬ 
tational  mesh  for  an  FE-ELLAM  discretization,  while  we  interpret  Mh  as  a 
point-distributed  grid  with  control  volumes  or  cells  around  each  of  the  interior 
vertices  Xi,i  =  1, . . . ,  ne-l,  so  that  0,  =  [xj_i/2,  Xi+1/2],  xi+i/2  =  (xi+xi+ 1)/2, 
and  | O.j |  =  Axi  =  xi+i/2  —  Xj-i/2-  In  the  following,  we  use  the  terms  mesh  and 
grid  interchangeably  to  refer  to  A4h  for  convenience. 

For  both  the  FE-ELLAM  and  FV-ELLAM,  the  trial  solution  is  represented 
using  the  standard  linear  Lagrangian  basis  functions 


ne 

C(x,  t )  Rj  c(x,  t)  =  YJ  Ci(t)i>i(x) 

i= 0 

where  C(x,t)  is  the  trial  solution  and  ^(x)  is  supported  on  [xj_i,Xj+i] 


(18) 


X  —  Xj-l 
Aa-i-i/2  ’ 


X  £  [Xi^Xi 


Xj+l—X 
A3-i+i/2  ’ 


X  £  [Xi,Xi+ 1] 


and  Axi+i/2  =  x%+\  —  x^  The  boundary  trial  functions  are 


(19) 


V’o  =  yr — a;e[a;L,a;i] 

AXi/2 

i  %  %ne  —  1  r  1  /on\ 

Vne  =  -7 - ,  X  £  [Xne-1,  XR\  (20) 

^xUe_  i/2 

where  O0  =  [^0^1/2]  and  =  [xne_i/2,  %ne\  at  the  boundaries  for  FV- 
ELLAM.  In  the  following,  we  write  C(x,t)  for  the  trial  solution  for  conve- 
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mence. 


At  time  level  tn+1,  the  test  function  Wj  is  aligned  in  a  regular  way  with  M.h. 
Since  wy  must  obey  eqn  (8),  this  will  not  be  the  case  in  general  for  t  <  tn+1. 
To  reinforce  this  observation,  we  use  the  notation  wf+1: 


<+10M"+I) 

<+10M”+1) 


i/ji,  for  i  —  0, . . . ,  ne  FE-ELLAM 

1,  X  G  [Xj_ i/2,  Xi+\/2[ 

FV-ELLAM 

0,  otherwise 


(21) 

(22) 


Similarly,  the  FV-ELLAM  boundary  test  functions  Wq+1  and  wn+1  are  indi¬ 
cator  functions  associated  with  fl0  and  QUe,  respectively. 


3. 1  Evaluation  of  integrals 


Given  a  choice  of  trial  and  test  functions,  the  next  step  is  to  approximate  the 
various  integrals  in  eqn  (17).  Since  we  use  a  backward  Euler  approximation  for 
the  dispersion  and  source  integrals  (second  and  last  terms  in  eqn  (17)),  these 
terms  can  be  evaluated  identically  to  similar  terms  that  appear  in  standard 
ELL  AM  formulations  for  linear  transport  problems  [8,  17]. 

3.1.1  mass  at  the  new  time  level 

The  first  term  in  eqn  (17)  accounts  for  the  mass  in  at  the  new  time  level. 
For  FE-ELLAM  we  have  the  standard  (nonlinear)  mass  integral  for  a  Galerkin 
finite  element  method 

J  M(Cn+1X+1  dx  for  %  =  0, . . . ,  ne  (23) 

a 

Similarly,  for  FV-ELLAM  we  have 
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(24) 


J  M  (C(x,  tn+1))w”+1(x,  tn+>)  dx—  f  M(C(x,  i"+1))dx 

Cl  Cli 

Eqns  (23)  and  (24)  can  be  approximated  numerically  in  a  straightforward  man¬ 
ner.  The  usual  approach  for  FV-ELLAM  discretizations  is  to  use  a  composite 
trapezoidal  rule,  which  is  exact  for  a  linear  problem  [17].  A  similar  approach 
can  be  used  here  for  both  the  FV-ELLAM  and  FE-ELLAM  discretizations, 
although  a  larger  number  of  subintervals  may  be  required  to  approximate 
the  nonlinear  mass  term  accurately.  Regardless  of  the  quadrature  used,  it  is 
also  important  for  mass  conservation  that  the  numerical  integration  strategies 
chosen  for  the  old  and  new  mass  integrals  be  consistent  [17]. 

3.1.2  mass  from  the  previous  time  level 

Tracking  along  characteristics  plays  a  major  role  in  the  remaining  integrals, 
which  account  for  boundary  contributions  and  the  mass  in  at  the  previ¬ 
ous  time  level.  As  a  result,  it  is  here  that  a  nonlinear  sorption  isotherm  can 
introduce  significant  complexity  over  conservative  transport  problems. 

Following  [15],  we  adopt  the  notation 

tn+! 

x*a(t)  =  Xa(t]x,tn+1)=x-  J  \a{C[x*(r),  t],x*(t),  r}  dr  (25) 

t 

t 

xa(t)  =  Xa(t;x,tn)  =x  +  jxa{C[xa(r),r],xa(r),r}dr  (26) 

tn 

for  tracking  along  adjoint  characteristics  given  by  eqn  (15).  For  tracking  along 
solution  characteristics  we  write 

tn+ 1 

x*(t)  —  X(t;  x,  tn+1)  =  x  —  j  \{C[x*(t),t],x*(t),t}  dr  (27) 

t 
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t 

x(t)  —  X(t]  x,  tn)  =  x  +  J  \{C[x(t),t],x(t),t}  dr  (28) 

tn 

Approximating  the  first  term  on  the  right  hand  side  of  eqn  (17)  numerically 
gives 


nq— 1 

M[C(x,  tn)]w?+1(x,  tn )  dx  «  £  M[Cn(xq)]w?+1(xq ,  tn)rq  (29) 

q= o 

where  we  define  the  set  of  integration  points  and  weights  {(xq,rq)},  q  = 
0, . . . ,  nq  —  1  in  physical  space  on  rather  than  on  a  reference  element.  nq  is 
the  total  number  of  integration  points.  We  employ  a  forward-tracking  approach 
where  the  integration  weights  and  points  are  defined  at  tn  to  ensure  accurate 
evaluation  of  the  mass  at  the  old  time  level.  The  mass  associated  with  each 
point  is  determined  by  evaluating  M[Cn{xq )]  at  the  interpolated  concentration 
values,  Cn{xq).  Each  integration  point  xq  is  tracked  forward  to  tn+1  by 

xq  =  (xq)a(tn+1)  =  Xa(tn+1;  xq,  tn)  (30) 

using  eqn  (26).  At  tn+1,  the  value  of  w™+1(xq,tn+1)  is  determined  using  eqn 
(21)  or  eqn  (22)  depending  upon  the  method  used. 

At  a  high  level,  approximating  the  mass  term  at  tn  simply  requires  iterating 
through  the  list  of  integration  points  xq  at  tn  and  tracking  them  forward  in 
time  to  tn+1  along  characteristics  defined  by  eqn  (26).  The  right  hand  side 
vector  components  of  eqn  (17)  for  test  functions  w™+1  that  are  nonzero  at  xq 
then  receive  a  corresponding  contribution  of  mass.  Unfortunately,  accurately 
tracking  the  characteristics  can  be  difficult  even  in  one  spatial  dimension.  As 
a  result,  a  number  of  strategies  have  been  developed  to  improve  the  tracking 
procedure’s  robustness.  In  particular,  FV-ELLAM  methods  introduce  an  ap- 
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proximate  test  function  W™+1  rather  than  use  the  actual  value  of  in  eqn 
(29).  Wi+l  is  roughly  a  smoothed  out  version  of  w™+1  with  wider  support  cov¬ 
ering  three  cells  in  one  spatial  dimension  and  is  intended  to  help  avoid  some 
of  the  difficulties  that  arise  in  distributing  mass  at  the  new  time  level  when 
tracking  is  inexact.  The  amount  of  smoothing  for  W™+1  is  given  in  terms  of 
the  parameter  NS  [17].  For  example,  in  the  interior  kF"+1  can  be  written 

x  <  Xi_  1/2  -  Aay_i 
Xi- 1/2  -  Aay_i  <  x  <  x{- l/2 
Xi- 1/2  <  X  <  Xi_  1/2  +  A  Xi 
Xi- 1/2  +  A  Xi  <X<  Xi+i/2  ~  A  Xi 
Xi+i/2  -  A  Xi  <  X  <  Xi+1/2 
Xi+1/2  A  X  <  Xi- 1-1/2  +  AXi+\ 

Xi+1/2  +  A%1  <  X 

(31) 

where  Ax  =  Aa;/NS.  FE-ELLAM  continues  to  use  w™+1  =  Note  that  for 
NS  =  2,  W™+ 1  is  also  piecewise-linear  chapeau  function.  However,  the  FV- 
ELLAM  and  FE-ELLAM  discretizations  differ  in  their  approximation  of  the 
new  mass  integral,  since  FV-ELLAM  continues  to  use  eqn  (24)  with  w^+1 
given  by  eqn  (22). 

We  use  the  composite  trapezoidal  rule  to  approximate  the  mass  integrals,  be¬ 
cause  it  has  superior  stability  properties  for  Eulerian-Lagrangian  discretiza¬ 
tions  [29]  and  allows  the  ability  to  adjust  the  number  of  subintervals  per  mesh 
cell  based  on  the  difficulty  of  a  given  problem  [17].  We  also  follow  the  common 
practice  for  FV-ELLAM  of  introducing  additional  integration  points  known  as 


Wi(x,t 


n+l  i 


0, 


Axi  I  1  _  xi— 1/2  x 

A:Ei_i+Axi  Ax;_i 

Axt+(x-Xj_1/2)Axi_i/A.Ti 

Axi-i+Axi 


=  < 


1, 

Axi+(xi+1/2-x)Axi+1/ Axj 


Axi+1+Axi 

Axz  f  -y  x~xi+ 1/2 

Axj+i+A xi  Axi+i 

0, 
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strategic  spatial  integration  points  (SSIPs)  and  strategic  temporal  integration 
points  (STIPs).  These  SSIPs  and  STIPs  are  obtained  by  backtracking  to  tn 
points  where  the  approximate  test  functions  W™+1  change  slope  and  so  im¬ 
prove  the  distribution  of  mass  at  the  new  time  level.  Defining  the  support  of 
W™+ 1  in  terms  of  NS  and  including  SSIPs  and  STIPs  has  proven  successful  for 
a  wide  range  of  applications  of  FV-ELLAM  [17,  31].  Usually,  the  number  of 
intervals  in  the  composite  trapezoidal  rule  is  also  determined  by  NS.  Since  the 
problems  here  are  nonlinear,  we  choose  the  number  of  subintervals  per  mesh 
cell  in  the  composite  trapezoidal  rule  based  on  accuracy  requirements  for  the 
nonlinear  mass  term  and  use  the  NS  parameter  only  to  define  the  support  of 
Wtn+1. 


Fig.  3.  £1*  (shaded)  Riemann  example  with  Cr  =  0 

In  the  following,  we  address  the  approximation  of  boundary  conditions  and 
detail  tracking  procedures  for  the  integration  points  {xq}.  Before  turning  to 
these  topics,  however,  the  behavior  of  solution  and  adjoint  characteristics  for 
zero  concentration  values  bears  some  additional  comment.  For  0  <  rif  <  1, 
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the  value  C  —  0  propagates  with  speed  0  even  with  non-zero  dispersion  [36] .  If 
we  consider  the  example  Riemann  problem  from  §2  with  Cr  =  0,  we  see  that 
Xa  =  cr.  Adjoint  characteristics  moving  with  speed  0  intersect  the  shock  and 
then  move  together  at  the  shock  speed  (see  Figure  3).  If  the  shock  intersects 
time  level  tn+l  at  x”+1,  then  Xa(tn;  x”+1,  tn+1)  is  non-unique.  Our  ELLAM 
approach  remains  well-defined  computationally  and  mathematically,  however. 
In  a  FV-ELLAM  context,  for  instance,  if  x™+1  £  O*  we  can  still  determine  the 
support  of  w"+1 

= [fe-1/2):,  fe+i/2):]  (32) 

from  the  points  (ay±i/2)*  =  Xa(tn;  (e"±i/2,  tn+1).  If,  say  x™+1  =  ay+ 1/2,  we  choose 
to  set  (xj+i /2)*  =  xi+i/2,  the  right-most  of  the  non-unique  values.  Moreover, 
since  M[C(xq,  tn)]  =  0  for  points  xq  at  tn  that  map  to  :r”+1,  the  integral  in  eqn 
(29)  is  independent  of  this  non-unique  choice.  Whenever  the  characteristics 
ahead  of  the  shock  carry  nonzero  mass  (Cr  >  0),  as  in  Figure  2,  they  do  not 
intersect,  and  Xa(tn]x™+1,tn+1)  is  unique. 

3.1.3  boundary  integrals 

The  third  integral  in  eqn  (17)  accounts  for  the  influence  of  the  physical  bound¬ 
ary.  Along  the  inflow  boundary,  eqn  (17)  contains  an  additional  term,  which 
for  T/  =  xl  is 


tn+1 

~  J  qb(xL,t)w(xLlt)  dt  (33) 

tn 

Eqn  (33)  results  in  a  contribution  to  the  right  hand  side  of  eqn  (17)  for  test 
functions  w™+1  that  intersect  the  inflow  boundary  over  [tn,tn+1].  Eqn  (33)  can 
be  approximated  much  as  eqn  (29)  where  the  numerical  quadrature  is  in  time 
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along  T j  [17].  Along  the  outflow  boundary  Y0  =  xR  the  zero-dispersive  flux 
condition  results  in  an  additional  term 

tU+l 

J  v(xR,t)C(xR,t)w(xR,t)  dt  (34) 

tn 

that  must  be  approximated.  We  simply  use  the  trapezoidal  rule 

tn+1 

J  v(xR,t)C(xR,t)dt  « 

tn 

A tn+1  r  1 

-j-  +U*h.«”+1)C”+1]  (35) 

since  this  requires  the  solution  from  time  levels  tn+1  and  tn  only  [31].  However, 
sub  time-stepping  along  To  can  be  used  for  greater  resolution  if  necessary. 

3.2  Tracking  techniques 

In  general,  a  numerical  integration  technique  is  required  to  solve  eqn  (14)  and 
determine  the  behavior  of  the  ELLAM  test  functions.  There  are  many  options 
for  linear  and  nonlinear  problems,  including  forward  and  backward  Euler  as 
well  as  explicit  Runge-Kutta  methods  [10].  To  track  from  point  (x°,r°)  to 
(x^r1)  with  a  backward  Euler  (BE)  approximation,  we  set 

xl  =  x°  +  Ar\a  [C{x\,  r1)]  (36) 

while  a  simple  second-order  explicit  Runge-Kutta  (RK2)  scheme  is  [10] 

%1a  =  x°  T  At Xa  [C (x°,  r0)] 
r1  =  r°  +  At 

xl  =  x°  +  {Aa[C'(^,T0)]  +  Aa[C'(^,T1)]}  (37) 

where  At  =  r1  —  r°  is  the  tracking  time  step. 
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The  ODE  integration  methods,  eqns  (36)  and  (37),  are  simple.  For  our  pur¬ 
poses,  the  main  challenge  in  applying  them  to  solve  eqn  (14)  is  the  evaluation 
of  A a[C(x,t)\.  Specifically,  when  integrating  over  the  time  interval  [tn,tn+1], 
the  value  of  the  approximate  solution  C(x,tn )  and  the  nonlinear  solver’s  cur¬ 
rent  guess  for  C(x,  tn+1)  will  be  known.  However,  determining  C(x,t)  for 
t  E  (tn,  tn+ 1 )  requires  interpolation  based  on  a  representation  over  0  x  [tn,  tn+1] 
for  the  approximate  solution.  There  are  a  variety  of  approaches  one  could  take 
to  obtain  such  a  representation  with  varying  degrees  of  complexity  depending 
on  the  assumptions  made  about  the  solutions  behavior.  If  one  wishes  to  use  a 
tracking  time  step  that  is  the  same  as  the  global  time  step,  then  eqn  (37)  or 
eqn  (36)  does  not  require  interpolation  at  intermediate  time  levels.  Of  course, 
this  suggests  that  the  Courant  number  allowed  may  then  be  limited  by  the 
accuracy  necessary  for  the  tracking  step.  Below,  we  denote  backward  Euler 
time  integration  with  At  =  A tn+1  as  BE-S  (backward  Euler,  single-step)  and 
label  the  combination  of  RK2  time  integration  with  At  =  A tn+1  as  RK2-S 
(second  order  Runge-Kutta,  single-step). 


We  also  consider  two  approaches  for  obtaining  intermediate  values  of  C  over 
(tn,tn+1]  to  allow  At  <  A tn+1.  The  first  is  to  use  bilinear  interpolation  (BL) 
in  space  and  time  based  on  C(x,tn)  and  the  nonlinear  solver’s  current  guess 
for  the  solution  at  the  new  time  level,  Cn+1,m(x,  tn+1).  The  second  approach  is 
to  use  a  front-tracking  algorithm  to  obtain  solution  estimates  at  intermediate 
time  levels,  tn,k  =  tn  +  kAr.  Specifically,  we  employ  a  front-tracking  method 
(FT)  based  on  Risebro  and  Tveito  [30]  and  implemented  in  Langseth  [26], 
which  computes  a  piecewise-constant  solution  to  the  homogeneous,  hyperbolic 
portion  of  eqn  (5), 
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(38) 


dM(C)  d(Cv) 
dt  dx 

We  then  use  the  solution  to  eqn  (38)  to  evaluate  A a[C(x,tn,k)\  in  eqn  (36) 
or  eqn  (37).  The  advantage  of  using  a  front-tracking  strategy  over  other  ap¬ 
proximate  solution  methods  is  that  it  can  accurately  locate  sharp  fronts  that 
cause  jumps  in  adjoint  characteristic  speeds  and  is  not  subject  to  a  Courant 
number  limitation  [26,  30].  We  note  that  the  front-tracking  method  itself  can 
not  be  applied  directly  to  the  original  nonlinear  transport  problem  due  to  the 
second-order  dispersion  term,  but  would  require  a  splitting  approach  [19,  21], 

The  initial  data  for  the  FT  front  tracking  could  in  principle  be  based  on 
either  C(x,tn )  or  on  Cn+1,m(x,tn+1).  Here,  we  use  C(x,tn )  and  track  forward 
to  tn+1.  To  evaluate  \a[C{x,  fn+1)]  we  use  the  tracking  solution  rather  than 
C'n+i’m(x,  tn+ 1 )  so  that  the  tracking  procedure  is  independent  of  the  solution  at 
the  new  time  level.  This  simplifies  the  Jacobian  calculation  and  nonlinear  solve 
dramatically,  but  eliminates  a  feedback  mechanism  present  in  the  tracking 
strategies  that  incorporate  Cn+1,m(x ,  tn+l).  Clearly,  the  front-tracking  strategy 
outlined  is  more  involved  than  bilinear  interpolation  or  relying  on  solution 
values  at  tn  and  tn+1  alone,  but  it  could  potentially  allow  significantly  larger 
Courant  numbers  and  so  fewer  nonlinear  solves  and  tracking  steps.  We  denote 
the  overall  tracking  procedure  RK2-FT  or  RK2-BL  when  either  the  FT  or  BL 
intermediate  solution  representation  is  combined  with  RK2  time  integration. 

We  mention  briefly  the  incorporation  of  boundary  data  in  the  tracking  proce¬ 
dures.  In  general,  one  can  simply  use  the  corresponding  value  from  the  trial 
solution  or  from  a  given  intermediate  solution  representation  when  evaluating 
Aa  at  boundary  locations.  For  Dirichlct  conditions,  the  concentration  at  the 
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inflow  boundary,  Cb,  can  be  used  when  evaluating  Aa.  This  may  also  be  a 
reasonable  approximation  for  total  flux  inflow  boundary  conditions  in  some 
cases,  such  as  advection-dominated  problems  where  an  inflow  concentration 
and  flow  rate  are  used  to  specify  qb. 


3.3  Nonlinear  and  linear  solvers 

The  FV-ELLAM  and  FE-ELLAM  approximations  presented  above  result  in 
a  discrete  nonlinear  system  at  each  time  level,  which  we  solve  using  Newton’s 
method.  This  nonlinear  solve  can  be  difficult  for  given  sets  of  physical  param¬ 
eters  and  auxiliary  data.  To  improve  robustness,  we  use  an  Armijo  line  search 
strategy  [22],  The  performance  of  Newton’s  method  can  also  be  improved  dra¬ 
matically  if  a  good  initial  guess  for  the  solution  is  available.  A  default  approach 
for  solution  of  transient  partial  differential  equations  with  implicit  time  dis¬ 
cretizations  is  to  use  the  solution  from  the  previous  time  step  as  the  initial 
guess  for  the  solution  at  the  new  time  level.  Obviously,  the  quality  of  this  guess 
degrades  as  the  size  of  the  time  step  increases  if  there  is  significant  transient 
behavior  in  the  problem.  To  improve  performance,  we  use  a  crude,  predicted 
value  based  on  linear  advection  only.  For  the  one- dimensional  problems  here, 
the  initial  guess  for  v  >  0  is, 

cf+1’°  =  q\  j  =  i-  LCrJ  (39) 

where  [xj  is  the  largest  integer  less  than  or  equal  to  x  and  Cr  is  the  target 
Courant  number  for  the  simulation  and  corresponds  to  the  maximum  char¬ 
acteristic  speed,  Cr  =  maxi=0(ie  |A(C'”)|Atn+1/Ax.  For  locations  that  track 
backwards  out  of  the  domain  (j  <  0  in  eqn  (39)),  we  set  Cf+1’°  to  the  cor- 
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responding  boundary  value.  The  intent  of  eqn  (39)  is  to  speed  convergence 
of  the  Newton  solve.  Its  impact  is  largest  for  the  RK2-S  and  BE-S  tracking 
strategies,  since  they  incorporate  values  of  Cn+1,m  in  the  tracking  of  numerical 
integration  points.  If  there  is  a  nonlinear  solver  failure,  the  time  step  is  halved 
and  the  Newton  solve  is  repeated  using  the  solution  from  the  last  unsuccessful 
iteration  as  the  initial  guess. 

Since  the  Freundlich  isotherm  is  not  differentiable  at  C  =  0  when  0  <  rif  <  1, 
we  evaluate  <p{C)  and  dp/dC  using  cubic  splines.  Whether  or  not  a  tracking 
strategy  depends  on  C(x,  tn+1)  has  a  significant  impact  on  the  complexity 
associated  with  calculating  analytical  Jacobians  for  the  Newton  solve.  We 
use  a  numerical  Jacobian  in  the  results  below  for  simplicity.  The  linear  solve 
needed  at  each  Newton  iteration  was  accomplished  using  a  banded  lower-upper 
decomposition  solver  from  LAPACK  [3]. 


4  Results 

We  next  present  a  series  of  numerical  experiments  to  evaluate  our  ELLAM 
approach  for  transport  problems  with  self-sharpening  fronts.  The  first  set  of 
simulations  examine  the  FV-ELLAM  and  FE-ELLAM  discretizations’  perfor¬ 
mance  with  the  RK2-S  and  BE-S  tracking  strategies  on  coarse  grids  and  the 
discretizations’  mass  conservation  properties.  We  then  investigate  the  meth¬ 
ods’  ability  to  resolve  fronts  accurately  as  solutions  become  more  steep.  Last, 
we  consider  the  performance  of  different  tracking  strategies  for  a  range  of  Cr. 

The  basic  test  problems  were  constant  injection  into  a  domain  originally  free 
of  contaminant  and  transport  of  a  contaminant  slug.  The  simulations  are  la- 
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beled  according  to  their  problem  (PA-PF),  which  ELLAM  discretization  was 
used  (FV  or  FE),  and  a  simulation  number.  In  all  cases,  the  porosity  and  ve¬ 
locity  were  6  =  0.4,  v  =  1.  The  final  simulation  time  was  t  =  0.5,  and  sorption 
was  modeled  by  a  Freundlich  isotherm  with  Kf  =  0.126  and  rif  =  0.7.  Zero 
dispersive  flux  was  applied  at  the  outflow  boundary.  For  the  injection  test 
problems,  the  inflow  boundary  condition  corresponded  to  a  constant  concen¬ 
tration  of  one.  Otherwise,  the  inflow  boundary  value  was  zero.  The  numerical 
quadrature  used  was  the  composite  trapezoidal  rule  with  six  subintervals.  The 
nonlinear  systems  in  all  simulations  were  solved  using  Newton’s  method  with 
a  numerical  Jacobian  and  an  I2  relative  residual  convergence  criterion.  The 
maximum  number  of  nonlinear  iterations  and  line  searches  allowed  was  twenty, 
and  the  nonlinear  solver  tolerance  was  ICE10  unless  stated  otherwise. 

4-1  Initial  results 

For  the  first  set  of  experiments,  we  considered  a  constant  injection  example 
(Problem  A)  and  transport  of  a  slug  initial  condition  (Problem  B)  given  by 

0,  0  <  x  <  0.15 

(x  —  0.15) /0.05,  0.15  <  x  <0.2 

C{x,t°)  =  *  l,  0.2  <  x  <0.3  (40) 

l-(x-  0.3) / 0.05,  0.3  <  x  <  0.35 

0,  0.35  <  x  <  1 

The  dispersion  coefficient  for  the  simulations  was  D  =  10”3.  The  other  relevant 
parameters  are  presented  in  Table  1,  where  Pe  is  the  mesh  Peclet  number,  and 
NT  is  the  number  of  composite  trapezoidal  rule  intervals  in  time  along  the 
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inflow  boundary.  We  note  that  Cr  is  the  target  Courant  number  for  a  given 
simulation  and  functions  as  an  upper  bound,  since  the  time  step  is  halved  if 
there  is  a  nonlinear  solver  failure. 

Table  1 

Run  parameters  for  initial  results 


Run 

Tracking 

Ax 

Pe 

Cr 

Run 

Tracking 

Ax 

Pe 

Cr 

PA.FV.l 

BE-S 

1/50 

20 

4.5 

PA.FV.6 

RK2-St 

1/100 

10 

4.5 

PA.FE.l 

BE-S 

1/50 

20 

4.5 

PA.FE.6 

RK2-St 

1/100 

10 

4.5 

PA.FV.2 

BE-S 

1/100 

10 

4.5 

PB.FV.l 

RK2-S 

1/50 

20 

4.5 

PA.FE.2 

BE-S 

1/100 

10 

4.5 

PB.FE.l 

RK2-S 

1/50 

20 

4.5 

PA.FV.3 

RK2-S* 

1/50 

20 

4.5 

PB.FV.2 

RK2-S 

1/100 

10 

4.5 

PA.FE.3 

RK2-S* 

1/50 

20 

4.5 

PB.FE.2 

RK2-S 

1/100 

10 

4.5 

PA.FV.4 

RK2-S* 

1/100 

10 

4.5 

PB.FV.3 

BE-S 

1/50 

20 

4.5 

PA.FE.4 

RK2-S* 

1/100 

10 

4.5 

PB.FE.3 

BE-S 

1/50 

20 

4.5 

PA.FV.5 

RK2-St 

1/50 

20 

4.5 

PB.FV.4 

BE-S 

1/100 

10 

4.5 

PA.FE.5 

RK2-St 

1/50 

20 

4.5 

PB.FE.4 

BE-S 

1/100 

10 

4.5 

NS  =  2;  Problem  A:  NT  =  16;  Problem  B:  NT  =  4 
*  BE-S  tracking  for  inflow  boundary 
f  RK2-S  tracking  using  A a(Cb)  at  inflow  boundary 


Figures  4-7  illustrate  the  performance  of  the  FV-ELLAM  and  FE-ELLAM  dis¬ 
cretizations  with  BE-S  tracking  for  Problem  A  and  RK2-S  tracking  for  Prob¬ 
lem  B.  The  corresponding  Li,  L2l  and  mass  balance  values  are  reported  in  Ta¬ 
bles  2  and  3.  The  dense  grid  solutions  were  obtained  with  a  mass  conservative 
finite  difference  discretization  on  grids  with  Ax  =  1/20000  and  Ax  =  1/50000, 
respectively.  For  these  simulations,  the  accuracy  of  the  FV-ELLAM  and  FE- 
ELLAM  discretizations  was  good.  The  FE-ELLAM  discretization  had  lower 
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Fig.  4.  ELLAM  solutions  Problem  A,  BE-S  tracking,  Ax  =  1/50,  t  =  0.5 


Fig.  5.  ELLAM  solutions  Problem  B,  RK2-S  tracking,  Ax  =  1/50,  t  =  0.5 
Li  and  L2  error,  but  had  some  overshoot  for  Run  PA.FE.l.  Otherwise,  both 
spatial  discretizations  were  able  to  represent  the  sharp  fronts  with  negligible 
over  or  undershoot.  The  mass  balance  results  were  good  for  both  methods  with 
the  FE-ELLAM  mass  error  at  the  level  of  the  nonlinear  solver  tolerance.  The 
FV-ELLAM  mass  balance  error  was  not  as  low  as  the  FE-ELLAM  error  due 
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Fig.  6.  ELLAM  solutions  Problem  A,  BE-S  tracking,  Ax  =  1/100,  t  =  0.5 


Fig.  7.  ELLAM  solutions  Problem  B,  RK2-S  tracking,  Ax  =  1/100,  t  =  0.5 
to  the  use  of  SSIPs  and  STIPs.  The  strategic  integration  points  increased  the 
FV-ELLAM  mass  balance  error  to  a  level  above  the  nonlinear  solver  residual 
because  the  numerical  integration  was  not  exact,  and  the  collection  of  integra¬ 
tion  points  changed  from  one  time  step  to  the  next  as  a  result  of  the  addition 
of  SSIPs  and  STIPs.  Global  mass  conservation  was  still  exact  up  to  the  non- 
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Table  2 


ELLAM  error  results  for  Problem  A 


Run 

Li 

L2 

rel.  mass  err 

PA.FV.l 

0.00514604 

0.0225308 

3.9826  x  10"6 

PA.FE.l 

0.00444420 

0.0186933 

— 

PA.FV.2 

0.00225567 

0.0110125 

1.8040  x  10"6 

PA.FE.2 

0.00177461 

0.00862980 

— 

PA.FV.3 

0.00440169 

0.0201363 

7.3856  x  10"6 

PA.FE.3 

0.00413805 

0.0171043 

— 

PA.FV.4 

0.00165866 

0.00839648 

7.9678  x  10-6 

PA.FE.4 

0.00130881 

0.00651466 

— 

PA.FV.5 

0.00442612 

0.0202831 

7.0136  x  10-6 

PA.FE.5 

0.00414199 

0.0171592 

— 

PA.FV.6 

0.00167139 

0.00842405 

8.3992  x  10-6 

PA.FE.6 

0.00130563 

0.00650926 

— 

—  less  than  10  10 


linear  solver  tolerance  for  each  step  from  tn  to  tn+1.  However,  the  strategic 
integration  points  were  defined  by  backtracking  from  the  new  time  level  lo¬ 
cations  where  the  approximate  test  functions  kF"+1  changed  slope.  For  this 
reason,  the  total  mass  in  the  domain  at  the  end  of  tn+1  defined  by  eqn  (24) 
was  in  some  cases  slightly  different  than  the  sum  in  eqn  (29)  for  the  step  from 
tn+1  to  tn+ 2  due  to  the  addition  or  subtraction  of  strategic  integration  points. 
This  effect  was  small,  however,  and  the  overall  solution  quality  was  improved 
by  the  use  of  the  STIP  and  SSIP  points.  We  also  note  that  the  FV-ELLAM 
discretization  achieved  second-order  convergence  in  space  for  Problem  B  for 
simulations  where  the  time  step  was  chosen  small  enough  to  eliminate  tempo- 
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Table  3 


ELLAM  error  results  for  Problem  B 


Run 

Li 

L2 

rel.  mass  err 

PB.FV.l 

0.00665642 

0.0208610 

7.6067  x  10“5 

PB.FE.l 

0.00520272 

0.0172889 

4.2842  x  10-10 

PB.FV.2 

0.00284313 

0.00917603 

1.9665  x  10"5 

PB.FE.2 

0.00228041 

0.00713305 

4.3568  x  10-10 

PB.FV.3 

0.0126638 

0.0287931 

1.1267  x  10“4 

PB.FE.3 

0.00735101 

0.0197214 

4.2842  x  10-10 

PB.FV.4 

0.00485904 

0.0130791 

4.4658  x  10-5 

PB.FE.4 

0.00336044 

0.00983110 

4.3568  x  10"10 

Fig.  8.  FV-ELLAM  solutions  Problem  A,  RK2-S  tracking  with  BE-S  tracking  at 
inflow,  Ax  =  1/50, 1/100,  t  =  0.5 

ral  truncation  error  [15]. 

In  Runs  PB.[FE,FV].3-4,  we  considered  BE-S  tracking  for  Problem  B.  The 
discretizations  still  resolved  the  fronts  with  negligible  overshoot  or  undershoot 
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and  produced  good  mass  balance  results,  but  the  Lx  errors  were  41%  and  47% 
higher  on  the  two  spatial  grids  for  FE-ELLAM  using  BE-S  tracking.  The  FV- 
ELLAM  L\  errors  with  BE-S  tracking  were  nearly  twice  as  large  as  the  RK2-S 
results.  Moreover,  the  PB.FE.3  and  PB.FE.4  simulations  experienced  one  and 
two  nonlinear  solver  failures,  respectively,  while  the  PB.FV.3  and  PB.FV.4 
runs  had  four  and  three  nonlinear  solver  failures,  respectively.  Reducing  the 
nonlinear  solver  tolerance  to  10“7  and  increasing  the  number  of  allowed  iter¬ 
ations  and  line  searches  to  100  did  not  change  the  number  of  nonlinear  solver 
failures  or  the  simulations’  basic  behavior. 

On  the  other  hand,  RK2-S  tracking  did  not  perform  as  well  for  Problem  A 
because  of  the  incompatibility  between  the  initial  data  and  the  non-zero,  to¬ 
tal  flux  inflow  boundary  condition.  The  Newton  solver  experienced  significant 
difficulties  when  Cn  was  used  as  the  initial  guess  for  the  solution  at  Cn+1. 
Even  with  the  initial  guess  from  eqn  (39),  both  FV-ELLAM  and  FE-ELLAM 
discretizations  exhibited  large  overshoot  at  early  times  when  the  trial  solution 
value  was  used  to  evaluate  Aa  at  the  inflow  boundary.  There  were  at  least  two 
approaches  to  improve  the  RK2-S  tracking  performance.  Given  the  advection- 
dominated  nature  of  problem  A,  we  considered  incorporating  a  boundary  con¬ 
centration,  Cb  =  1,  into  the  RK2-S  tracking  by  using  \a(Cb)  for  points  at 
the  inflow  boundary,  xl ■  Another  solution  was  simply  to  use  RK2-S  tracking 
for  interior  points  and  BE-S  tracking  at  the  inflow  boundary.  The  simulations 
with  the  combined  tracking  for  Problem  A  are  PA.[FE,FV].3-4.  Simulations 
using  A a(Cb)  for  inflow  boundary  values  are  labeled  PA.[FE,FV].5-6.  The  ap¬ 
proaches  performed  similarly  and  eliminated  the  nonlinear  solver  difficulties. 
As  the  results  in  Table  2  and  Figure  8  indicate,  both  strategies  led  to  improved 
accuracy  over  BE-S  tracking  alone  for  the  FE-ELLAM  and  FV-ELLAM  spa- 
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tial  discretizations. 


4-2  Front  resolution  examples 


In  the  runs  for  Problem  A  and  Problem  B,  both  the  FV-ELLAM  and  FE- 
ELLAM  discretizations  were  able  to  resolve  the  contaminant  fronts  adequately 
for  a  mesh  width  of  Ax  =  1/50,  where  the  solution  fronts  spanned  approx¬ 
imately  three  intervals.  Given  the  use  of  a  continuous,  piecewise-linear  rep¬ 
resentation  for  the  trial  solution,  this  is  essentially  the  minimum  we  could 
expect  and  agrees  with  previous  results  for  ELLAM  and  MMOC  discretiza¬ 
tions  of  linear  transport  problems  [33].  For  conditions  where  sharper  fronts 
arise,  say  with  less  physical  dispersion  or  lower  Freundlich  exponents,  one  can 
expect  that  the  ELLAM  discretizations  would  require  finer  grids  to  resolve 
the  solution  monotonically. 

To  be  more  specific,  we  looked  at  the  discretizations’  performance  as  the  front 
widths  in  our  test  problems  decreased.  A  straightforward  way  to  do  this  was 
to  look  at  traveling  wave  solutions  to  eqns  (5)-(6)  with  boundary  data  corre¬ 
sponding  to  Problem  A.  That  is,  we  sought  a  solution  C(rj)  to  eqns  (5)-(6) 
over  the  real  axis  with  boundary  data  C(— 00)  =  C*  and  C( 00)  =  C*  where 
rj  —  x  —  at  and  a  is  the  wave  speed  [37].  For  equilibrium  sorption  with  a 
Freundlich  isotherm  and  C*  =  0,  van  Duijn  and  Knabner  [38]  provide  a  closed 
form  expression  for  the  traveling  wave  solution 


C(v)  =  c* 

0, 


1 1  —  exp 
otherwise 


(1  -  nf) 


for  rj  <  0 


where  the  wave  speed  is 


(41) 
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c*  -  a 


a  —  v 


ip(C*)  -  <p(C *)  +  c*  -  c* 


(42) 


The  solution  front  width  can  be  determined  by  inverting  eqn  (41)  for  rj 


inq-tc/c-)1^  for0<c<c. 


=  D 


(1-v)  L?) 

In  (1  -  (C'/C'*)1“n/) 


(43) 


(1  -  nf)  (■ v  -  a) 

and  calculating  the  distance  from  the  point  C  =  0  to  the  point  where  C  = 
C*(  1  —  e),  since  the  solution  to  the  problem  is  non- increasing.  The  front  is 
located  at  rj  =  0,  so  we  set  the  front  width  to  be  di_e  =  \rj[C*(l  —  e)]|.  The 
front  width  for  a  fixed  e,  do.99  f°r  example,  scales  linearly  with  the  disper¬ 
sion  D,  while  the  front  scales  with  \J~D  at  a  fixed  time  for  the  solution  to 
the  corresponding  linear  advection-dispersion  problem  [33].  In  particular,  eqn 
(43)  suggests  that  if  our  ELLAM  approximation  needs  roughly  three  inter¬ 
vals  to  resolve  a  front  without  overshoot  or  undershoot,  then  the  minimum 
discretization  required  will  be  given  by 


^0.99/3 


(44) 


Table  4  gives  run  parameters  for  ELLAM  calculations  with  D  =  2.5  x  10-4, 
Ax  =  1/100, 1/200,  and  an  initial  condition  given  by  eqn  (41).  Aside  from  D, 
the  physical  parameters  were  identical  to  those  for  Problem  A.  The  front  width 
was  do. 99  =  1-4498  x  10“2,  so  that  do.99/3  =  4.8327  x  10^3.  Figure  9  shows  the 
solutions  for  Ax  =  1/100  and  Figure  10  shows  the  results  for  Ax  =  1/200. 
The  corresponding  error  values  are  reported  in  Table  5. 

For  Problem  C,  we  performed  simulations  with  FV-ELLAM  using  two  different 
values  of  NS.  NS  =  2  represents  a  conservative  choice,  since  the  approximate 
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Table  4 


Run  parameters  for  Problem  C 


Run 

Tracking 

Ax 

D 

Pe 

Cr 

NS 

PC.FV.l 

RK2-S* 

1/100 

2.5  x  10-4 

40 

4.5 

2 

PC.FE.l 

RK2-S* 

1/100 

2.5  x  10"4 

40 

4.5 

- 

PC.FV.2 

RK2-S* 

1/200 

2.5  x  10"4 

20 

4.5 

2 

PC.FE.2 

RK2-S* 

1/200 

2.5  x  10-4 

20 

4.5 

- 

PC.FV.3 

RK2-S* 

1/100 

2.5  x  10"4 

40 

4.5 

16 

PC.FV.4 

RK2-S* 

1/200 

2.5  x  10-4 

20 

4.5 

16 

NT  =  16 


*  BE-S  tracking  for  inflow  boundary 
Table  5 

ELLAM  error  results  for  Problem  C 


Run 

Li 

L-2 

rel.  mass  err 

PC.FV.l 

0.00307715 

0.0231885 

3.9674  x  10"5 

PC.FE.l 

0.00321725 

0.0218095 

1.5572  x  10-10 

PC.FV.2 

0.000984260 

0.00967998 

9.0764  x  10"6 

PC.FE.2 

0.000921269 

0.00836204 

1.7097  x  10"9 

PC.FV.3 

0.00443214 

0.0305197 

1.7577  x  10~4 

PC.FV.4 

0.000906421 

0.00828032 

2.5945  x  10-5 

test  function  bE/1"1"1  is  a  piecewise-linear  chapeau  function.  The  wider  support 
for  NS  =  2  has  the  effect  of  spreading  mass  over  a  larger  area  in  the  approx¬ 
imation  of  the  old  mass  integral,  eqn  (29).  Choosing  NS  =  16  narrows  the 
support  for  W™+1  significantly  and  leads  to  a  less  diffusive  approximation.  As 
Figure  9  shows,  the  difference  in  the  solutions  with  NS  =  2  and  NS  =  16 
for  Aa;  =  1/100  was  dramatic.  However,  for  the  finer  grid  with  Ax  ~  d0  .99/3 
both  approximations  resolved  the  solution  front  accurately  without  spurious 
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Fig.  9.  ELLAM  results  for  Problem  C,  Ax  =  1/100,  t  =  0.5 


oscillations. 


The  FV-ELLAM  approximate  test  functions  W™+1  are  equivalent  to  FE- 
ELLAM  test  functions  when  NS=2.  However,  W™+1  is  only  used  in  the  FV- 
ELLAM  approximation  of  mass  at  tn:  so  the  discretizations  are  not  identical. 
More  precisely,  the  FE-ELLAM  approximation  to  the  mass  at  tn+l ,  eqn  (23), 
corresponds  to  a  more  distributed  mass  matrix  than  the  FV-ELLAM  approx¬ 
imation,  eqn  (24),  and  so  contains  less  numerical  diffusion  [17,  25,  31].  As  a 
result,  the  the  FE-ELLAM  approximation  produced  visible  overshoot  for  the 
coarser  grid  but  still  resolved  the  front  accurately  for  Ax  ~  do. 99/3. 

We  also  performed  simulations  for  the  slug  initial  condition  with  Ax  =  1/100, 
1/200,  and  D  =  2.5  x  10~4.  While  not  a  traveling  wave  solution,  the  front 
behavior  was  similar  to  Problem  C.  Table  6  gives  run  parameters  for  Problem 
D  and  Table  7  gives  the  error  values  for  the  ELLAM  calculations.  Figures 
11  and  12  show  the  results  for  the  two  grids.  The  dense  grid  solution  was 
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Fig.  10.  ELLAM  results  for  Problem  C,  Ax  =  1/200,  t  =  0.5 

again  obtained  with  a  mass  conservative  finite  difference  discretization  on  a 
grid  with  Ax  =  1/50000.  The  relative  performance  of  the  discretizations  was 
similar  to  Problem  C.  All  the  discretizations  produced  accurate  solutions  for 
Ax  ps  Axmm ,  but  the  FE-ELLAM  and  the  FV-ELLAM  discretization  with 
NS  =  16  exhibited  overshoot  on  the  Ax  =  1/100  grid. 

Table  6 

Run  parameters  for  Problem  D 


Run 

Tracking 

Ax 

D 

Pe 

Cr 

NS 

PD.FV.l 

RK2-S 

1/100 

2.5  x  10~4 

40 

4.5 

2 

PD.FE.l 

RK2-S 

1/100 

2.5  x  10-4 

40 

4.5 

- 

PD.FV.2 

RK2-S 

1/200 

2.5  x  10“4 

20 

4.5 

2 

PD.FE.2 

RK2-S 

1/200 

2.5  x  10”4 

20 

4.5 

- 

PD.FV.3 

RK2-S 

1/100 

2.5  x  10"4 

40 

4.5 

16 

PD.FV.4 

RK2-S 

1/200 

2.5  x  10"4 

20 

4.5 

16 

NT  =  4 
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Table  7 


ELLAM  error  results  for  Problem  D 


Run 

L  i 

L2 

rel.  mass  err 

PD.FV.l 

0.00375488 

0.0222302 

2.9905  x  10“5 

PD.FE.l 

0.00373496 

0.0211762 

- 

PD.FV.2 

0.00127395 

0.00884059 

1.6434  x  10”5 

PD.FE.2 

0.00115919 

0.00777404 

8.7131  x  10"10 

PD.FV.3 

0.00398739 

0.0218459 

9.6868  x  10“5 

PD.FV.4 

0.00114844 

0.00777340 

3.5685  x  10“5 

-  less  than  10  10 


Fig.  11.  ELLAM  results  for  Problem  D,  Ax  =  1/100,  t  =  0.5 
4-3  Tracking  comparison 

The  simulations  for  Problems  A-D  were  performed  using  the  RK2-S  and  BE-S 
tracking  strategies,  which  required  Ar  =  A tn+1.  While  we  were  able  to  obtain 
good  results  for  Cr  =  4.5,  relying  solely  on  characteristic  speeds  from  tn  and 
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Fig.  12.  ELLAM  results  for  Problem  D,  Ax  =  1/200,  t  =  0.5 

tn+1  is  less  tenable  as  Cr  increases  significantly.  We  next  compared  the  RK2-S 
tracking  to  the  RK2-BL  and  RK2-FT  strategies  with  At  =  A tn+1/nst,  where 
the  number  of  tracking  substeps,  nst ,  was  greater  than  one.  First,  we  compared 
the  RK2-FT  and  RK2-BL  strategies  for  moderate  Cr  where  the  RK2-S  strat¬ 
egy  performed  well.  The  physical  parameters  for  the  first  set  of  simulations 
(Problem  E)  were  identical  to  Problem  B.  We  then  used  the  physical  param¬ 
eters  and  initial  condition  from  Problem  D  to  compare  the  tracking  strategies 
for  significantly  larger  Cr  (Problem  F). 


Tables  8  and  9  contain  the  relevant  parameters  for  the  simulations  of  Problem 
E  and  Problem  F,  respectively.  The  main  parameters  controlling  the  perfor¬ 
mance  of  the  RK2-FT  method  are  the  discretization  parameter  5fr  and  the 
wave  speed  cutoff,  c/r.  5fr  determines  the  number  of  fronts  used  to  represent 
the  solution  initially  and  the  number  of  shock  waves  used  to  approximate  a 
rarefaction.  cjr  is  used  along  with  Sfr  to  provide  a  minimum  wave  strength 
required  for  a  front  to  be  incorporated  in  the  solution  approximation.  Specih- 
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cally,  a  front  is  tracked  if  the  magnitude  of  the  jump  across  the  front  is  greater 
than  Cfr5fr  [26]- 
Table  8 

Run  parameters  for  Problems  E 


Run 

Tracking 

Ax 

D 

Pe 

Cr 

nst 

5fr 

Gfr 

PE.FV.l 

RK2-BL 

1/50 

1.0  x 

10-3 

20 

4.5 

2 

- 

- 

PE.FE.l 

RK2-BL 

1/50 

1.0  x 

10"3 

20 

4.5 

2 

- 

- 

PE.FV.2 

RK2-FT 

1/50 

1.0  x 

10-3 

20 

4.5 

2 

1/50 

0.01 

PE.FE.2 

RK2-FT 

1/50 

1.0  x 

10~3 

20 

4.5 

2 

1/50 

0.01 

PE.FV.3 

RK2-BL 

1/100 

1.0  x 

10-3 

10 

4.5 

2 

- 

- 

PE.FE.3 

RK2-BL 

1/100 

1.0  x 

10"3 

10 

4.5 

2 

- 

- 

PE.FV.4 

RK2-FT 

1/100 

1.0  x 

10-3 

10 

4.5 

2 

1/100 

0.01 

PE.FE.4 

RK2-FT 

1/100 

1.0  x 

10-3 

10 

4.5 

2 

1/100 

0.01 

PE.FV.5 

RK2-FT 

1/100 

1.0  x 

10"3 

10 

4.5 

2 

1/50 

0.01 

PE.FE.5 

RK2-FT 

1/100 

1.0  x 

10-3 

10 

4.5 

2 

1/50 

0.01 

PE.FV.6 

RK2-FT 

1/100 

1.0  x 

10"3 

10 

4.5 

2 

1/200 

0.01 

PE.FE.6 

RK2-FT 

1/100 

1.0  x 

10~3 

10 

4.5 

2 

1/200 

0.01 

NS  =  2,  NT  =  4 


Figure  13  shows  FV-ELLAM  solutions  from  the  RK2-BL  and  RK2-FT  track¬ 
ing  strategies  for  Problem  E  with  Ax  =  1/50  and  Cr  =  4.5.  For  comparison, 
the  solution  using  RK2-S  tracking,  Run  PB.FV.l,  is  included  as  well.  Table 
10  includes  the  L1;  L2,  and  mass  balance  errors  for  the  RK2-BL  and  RK2-FT 
tracking  strategies.  The  results  for  the  RK2-BL  tracking  were  much  poorer 
than  either  the  RK2-S  or  the  RK2-FT  methods,  producing  large  overshoot  on 
both  the  Ax  =  1/50  and  Ax  =  1/100  grids.  The  RK2-FT  solutions  were  less 
accurate  than  the  RK2-S  results  for  the  same  Cr  for  both  the  FE-ELLAM 
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Table  9 


Run  parameters  and  results  for  Problems  F 


Run 

Tracking 

Ax 

D 

Pe 

Cr 

nst 

5fr 

cfr 

PF.FE.l 

RK2-FT 

1/200 

2.5  x  10~4 

20 

16.5 

2 

1/200 

0.01 

PF.FV.l 

RK2-FT 

1/200 

2.5  x  10~4 

20 

16.5 

2 

1/200 

0.01 

PF.FE.2 

RK2-S 

1/200 

2.5  x  10"4 

20 

16.5 

1 

- 

- 

PF.FV.2 

RK2-S 

1/200 

2.5  x  10”4 

20 

16.5 

1 

- 

- 

PF.FE.3f 

RK2-S 

1/200 

2.5  x  10”4 

20 

16.5 

1 

- 

- 

PF.FV.3t 

RK2-S 

1/200 

2.5  x  10”4 

20 

16.5 

1 

- 

- 

PF.FE.4 

RK2-FT 

1/200 

2.5  x  10”4 

20 

16.5 

2 

1/100 

0.01 

PF.FV.4 

RK2-FT 

1/200 

2.5  x  10-4 

20 

16.5 

2 

1/100 

0.01 

NT  =  4 


f  nonlinear  solver  tolerance  10  6,  max  100  line  searches,  100  nonlinear  iterations 

and  FV-ELLAM  spatial  discretizations.  While  fronts  were  similarly  resolved, 
the  L\  error  for  the  RK2-FT  Runs  was  35%  higher  for  FV-ELLAM  and  19% 
higher  for  FE-ELLAM  on  the  Ax  =  1/50  grid. 

The  impact  of  5fr  on  the  accuracy  of  the  RK2-FT  tracking  can  be  seen  in 
Runs  PE.[FV,FE].4-PE.[FV,FE].6.  RK2-FT  tracking  with  5fr  =  1/100  and 
Sfr  =  1/200  produced  Li  and  L-2  errors  that  were  comparable  to  the  errors  from 
the  corresponding  RK2-S  tracking  solutions.  Increasing  5fr  to  2 Ax  produced 
less  accurate  results,  as  can  be  seen  in  Figure  14  and  Table  10. 

Since  the  RK2-BL  strategy  performed  poorly  for  Problem  E,  we  considered 
only  the  RK2-S  and  RK2-FT  approaches  for  the  problem  with  less  physi¬ 
cal  dispersion.  Figure  15  shows  the  solutions  from  the  RK2-S  and  RK2-FT 
strategies  with  FV-ELLAM  and  FE-ELLAM  spatial  discretizations  for  a  tar- 
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Table  10 


Run  results  for  Problems  E 


Run 

Li 

l2 

rel.  mass  err 

PE.FV.l 

0.0155545 

0.0383824 

4.1861  x  10“5 

PE.FE.l 

0.0160688 

0.0396553 

- 

PE.FV.2 

0.00899326 

0.0240897 

6.6298  x  10"5 

PE.FE.2 

0.00620536 

0.0188605 

- 

PE.FV.3 

0.00674198 

0.0184749 

2.4097  x  10“5 

PE.FE.3 

0.00711775 

0.0191097 

4.3568  x  10~10 

PE.FV.4 

0.00359544 

0.0112651 

2.7236  x  10-5 

PE.FE.4 

0.00235988 

0.00807344 

8.7136  x  10"10 

PE.FV.5 

0.00802655 

0.0153431 

3.1588  x  10-5 

PE.FE.5 

0.00698441 

0.0132661 

4.3568  x  10”10 

PE.FV.6 

0.00334196 

0.0111522 

1.7654  x  10-5 

PE.FE.6 

0.00221380 

0.00813326 

- 

-  less  than  10  10 


get  Cr  =  16.5,  while  Table  11  gives  the  L1;  L2,  and  mass  balance  errors  as 
well  as  the  total  number  of  time  steps  taken  and  nonlinear  solver  iterations  for 
the  simulations.  The  spatial  discretizations  performed  similarly.  Results  with 
both  tracking  approaches  captured  the  sharp  front  well  with  negligible  mass 
balance  error. 

The  Li  and  L>  errors  for  the  RK2-S  and  RK2-FT  simulations  with  Sfr  —  Ax 
were  similar.  The  FE-ELLAM  errors  with  RK2-S  tracking  were  lower  than 
those  for  FE-ELLAM  with  RK2-FT  tracking,  while  the  opposite  was  true  for 
the  FV-ELLAM  discretization.  The  RK2-FT  results  with  Sfr  —  2 Ax  were 
much  poorer.  Although  the  accuracy  in  Runs  PF.[FE,FV].2  was  comparable 
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Table  11 


Run  results  for  Problems  F 


Run 

Li 

rel.  mass  err 

time 

steps 

nl  its 

PF.FV.l 

0.00217702 

0.0129708 

3.2722  x  10"6 

5 

27 

PF.FE.l 

0.00196901 

0.0122595 

1.3070  x  10~9 

5 

28 

PF.FV.2 

0.00240612 

0.0106530 

1.6912  x  10"4 

54 

462 

PF.FE.2 

0.00134290 

0.00757244 

1.7426  x  10-9 

46 

425 

PF.FV.3 

0.00189700 

0.0102873 

1.9320  x  10~4 

36 

362 

PF.FE.3 

0.00134111 

0.00757331 

6.8370  x  10-5 

53 

433 

PF.FV.4 

0.00594223 

0.0162554 

9.1099  x  10"6 

5 

27 

PF.FE.4 

0.00592168 

0.0161490 

- 

5 

30 

-  less  than  10  10 


Fig.  13.  FV-ELLAM  results  for  Problem  E,  RK2-S,  RK2-BL,  and  RK2-FT  tracking, 
Ax  =  1/50,  t  =  0.5 
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Fig.  14.  FE-ELLAM  results  for  Problem  E,  RK2-FT  tracking  8fr  =  1/50,1/100, 
Ax  =  1/100,  t  =  0.5 


to  that  in  Runs  PF.[FE,FV].l,  Cr  =  16.5  was  too  large  for  the  RK2-S  tracking. 
Both  Runs  PF.FV.2  and  PF.FE.2  experienced  repeated  nonlinear  solver  fail¬ 
ures  due  to  the  larger  target  Cr.  The  FV-ELLAM  simulation  had  16  failures 
and  took  54  time  steps,  while  the  FE-ELLAM  calculation  had  14  failures  and 
46  time  steps  overall.  The  total  number  of  nonlinear  iterations  was  an  order  of 
magnitude  higher  as  well.  The  failure  modes  for  the  Newton  solver  consisted 
of  both  stagnation  in  the  Armijo  line  search  and  exhausting  the  allowed  non¬ 
linear  iterations.  Since  the  relative  residual  tolerance  of  10“10  was  rather  tight, 
Runs  PF.FV.3  and  PF.FE.3  were  performed  with  a  tolerance  of  10-6  and  a 
maximum  of  100  nonlinear  iterations  and  line  searches.  The  reduced  tolerance 
decreased  the  number  of  failures  for  FV-ELLAM  to  13  but  increased  the  num¬ 
ber  of  failures  for  the  FE-ELLAM  simulation  to  15.  The  number  of  time  steps 
and  nonlinear  iterations  required  were  still  an  order  of  magnitude  higher  for 
the  RK2-S  tracking  simulations. 
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Fig.  15.  ELLAM  results  for  Problem  F,  RK2-S  and  RK2-FT  tracking,  Ax  =  1/200, 
t  =  0.5 


5  Discussion 


The  test  problems  we  considered  were  challenging  because  they  involved  self- 
sharpening  fronts  with  small  physical  dispersion  as  well  as  rarefaction  tails. 
While  there  were  some  differences,  the  performance  and  accuracy  of  the  FV- 
ELLAM  and  FE-ELLAM  discretizations  were  comparable  in  the  numerical 
experiments.  Both  discretizations  resolved  the  fronts  accurately  as  long  as  the 
mesh  width  was  sufficient  to  allow  three  elements  on  a  front.  FE-ELLAM  gen¬ 
erally  produced  lower  Li  and  L2  error  but  was  more  prone  to  overshoot  or 
undershoot  on  coarse  grids  than  FV-ELLAM  with  NS=2.  With  NS=2,  our 
FV-ELLAM  approach  had  more  numerical  diffusion  than  FE-ELLAM  due 
to  the  differences  in  their  approximations  for  the  mass  at  tn+1  [17,  25,  31]. 
Higher  values  of  NS  reduced  the  numerical  diffusion  for  FV-ELLAM  and  im¬ 
proved  accuracy  for  sufficiently  resolved  calculations,  but  naturally  increased 
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the  likelihood  of  spurious  oscillations  on  coarse  grids.  In  general,  there  was 
the  familiar  tradeoff  between  accuracy  in  smooth  regions  and  the  ability  to 
represent  fronts  monotonically  [27].  Steeper  physical  solutions  required  hirer 
discretizations  or  an  approximation  with  greater  numerical  diffusion  to  create 
fronts  with  a  minimum  width  for  a  given  grid. 

The  global  mass  balance  errors  for  the  FE-ELLAM  and  FV-ELLAM  discretiza¬ 
tions  were  good.  The  FE-ELLAM  mass  balance  error  was  at  the  nonlinear 
solver  tolerance,  while  the  FV-ELLAM  errors  were  below  2  x  10~4  for  all  the 
simulations.  The  FV-ELLAM  mass  balance  error  was  not  as  low  as  the  FE- 
ELLAM  error  due  to  the  use  of  SSIPs  and  STIPs  [15].  On  the  other  hand,  the 
SSIPs  and  STIPs  improved  the  overall  performance  of  the  FV-ELLAM  dis¬ 
cretization,  and  the  FV-ELLAM  relative  mass  balance  errors  remained  small. 
The  impact  of  the  strategic  integration  points  on  global  mass  balance  can  be 
reduced  by  increasing  the  accuracy  of  the  quadrature  through,  for  instance, 
adding  more  intervals  in  the  composite  trapezoidal  rule  formulas. 

We  explored  two  time  discretizations,  BE  and  RK2,  for  tracking  numerical 
integration  points  forward  from  tn  to  tn+ 1 .  The  first-order  BE  method  was  not 
as  accurate  as  the  second-order  RK2  approximation  when  using  characteristic 
information  from  tn  and  tn+l  alone.  On  the  other  hand,  it  was  useful  for 
tracking  from  the  inflow  boundary  when  qb  >  0.  One  of  the  advantages  of 
ELLAM  is  that  the  tracked  integration  points  are  independent  of  one  another, 
so  it  was  then  simple  to  combine  different  tracking  procedures  within  the  same 
simulation. 

Since  tracking  required  information  about  the  solution  to  evaluate  the  ad¬ 
joint  characteristics,  we  also  investigated  different  approaches  for  obtaining 
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a  solution  representation  at  times  t  >  tn.  The  RK2-S  approach,  which  was 
simple  and  relied  solely  on  Cn  and  Cn+ 1  ,m  worked  well  for  Cr  up  to  4.5  in  the 
numerical  experiments  when  combined  with  the  initial  guess  from  eqn  (39). 
In  general,  an  upper  bound  on  the  time  step  for  the  RK2-S  approach  de¬ 
pends  on  the  difficulty  of  the  problem  and  the  spatial  discretization,  including 
the  amount  of  numerical  diffusion.  For  time  steps  that  were  too  large,  poor 
performance  of  the  RK2-S  tracking  algorithm  manifested  itself  through  non¬ 
linear  solver  failures  for  the  target  time  step,  which  in  turn  led  to  time  step 
reductions  in  order  to  obtain  convergence.  On  the  other  hand,  the  RK2-FT 
tracking  strategy,  which  combined  RK2  time  integration  and  a  front-tracking 
method  to  provide  intermediate  solution  values,  allowed  the  FE-ELLAM  and 
FV-ELLAM  methods  to  obtain  good  solutions  for  very  large  time  steps. 

Accuracy  of  the  intermediate  representation  naturally  affected  the  overall  so¬ 
lution  accuracy,  since  the  tracking  dictated  the  mass  distribution  from  the  pre¬ 
vious  time  step.  Using  bilinear  interpolation  to  obtain  intermediate  solution 
representations  (RK2-BL)  performed  poorly.  For  smaller  time  steps  where  the 
RK2-S  tracking  succeeded,  the  RK2-FT  solutions  generally  had  higher  L\  and 
L2  error  values  than  those  from  the  RK2-S  tracking  algorithm.  Simulations 
with  the  two  tracking  strategies  resolved  fronts  similarly.  Rather,  the  increased 
error  for  the  RK2-FT  approximation  was  largely  around  rarefactions.  Increas¬ 
ing  Sfr  improved  the  FT  algorithm’s  approximation  of  rarefactions  and  in  turn 
reduced  the  solution  error. 

In  our  numerical  experiments,  both  the  FV-ELLAM  and  FE-ELLAM  dis¬ 
cretizations  resolved  the  self-sharpening  fronts  monotonically  as  long  as  the 
mesh  width  was  sufficient  to  allow  three  intervals  on  a  front.  This  places  an 
increasing  computational  burden  on  the  discretizations  as  fronts  in  a  prob- 
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lem  steepen  clue,  for  example,  to  reduced  physical  dispersion.  On  the  other 
hand,  increased  numerical  diffusion  can  always  be  introduced  using  techniques 
such  as  mass  lumping  [31,  34]  in  order  to  widen  solution  fronts  sufficiently  to 
meet  the  requirements  on  a  given  mesh.  Another  alternative  to  improve  per¬ 
formance  is  to  employ  adaptive  refinement  around  the  front  using  patch-based 
local  refinement  [39]  or  a  moving  mesh  approach  [43,  44]  in  order  to  provide 
the  minimal  resolution  necessary  around  fronts  with  coarser  discretizations  in 
areas  of  smoother  solution  behavior. 

We  have  not  addressed  the  issue  of  computational  efficiency  in  this  work. 
For  multidimensional  nonlinear  problems,  we  expect  the  expense  and  suc¬ 
cess  of  our  approach  to  be  dictated  largely  by  the  tracking  procedure  for 
numerical  integration  points.  Here,  we  have  identified  promising  strategies  for 
one-dimensional  problems.  The  RK2-S  tracking  is  simple  to  implement  and 
inexpensive  on  a  per  iteration  basis.  The  RK2-FT  is  more  involved  since  it 
requires  the  initialization  and  solution  of  a  front-tracking  problem  every  time 
step.  However,  the  RK2-FT  allowed  much  larger  time  steps  than  the  RK2-S 
tracking  for  the  problems  considered.  Another  aspect  of  the  RK2-FT  tracking 
procedure,  as  implemented  here,  was  that  it  was  independent  of  the  solution  at 
the  new  time  level.  While  this  was  not  necessary,  it  did  simplify  the  Jacobian 
calculation  and  behavior  of  the  nonlinear  systems  dramatically,  which  could 
offer  a  significant  advantage  for  multidimensional  problems. 

While  our  initial  results  are  promising  given  the  straightforward  nature  of  the 
formulation,  more  work  is  required  to  evaluate  its  performance  for  systems 
which  present  significantly  more  complicated  characteristic  behavior  and  for 
multidimensional  problems.  For  problems  in  two  and  three  spatial  dimensions, 
much  of  the  current  approach  can  be  extended  naturally  following  previous 
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ELLAM  work  [8,  18].  The  main  challenge  will  be  to  develop  effective  tracking 
procedures  given  the  added  nonlinearity  of  the  adjoint  characteristics. 


6  Conclusions 


In  this  work,  we  formulated  FV-ELLAM  and  FE-ELLAM  approximations  for 
advective-dispersive  transport  with  nonlinear  equilibrium  sorption.  We  per¬ 
formed  a  series  of  numerical  experiments  to  evaluate  the  discretizations’  front 
resolution,  mass  balance  properties,  and  ability  to  take  large  time  steps  while 
maintaining  accuracy.  Based  on  our  work,  we  draw  the  following  conclusions: 

•  Our  ELLAM  formulation  produces  solutions  that  maintain  the  mass-conservation 
properties  for  which  ELLAMs  are  known.  Both  FV-ELLAM  and  FE-ELLAM 
discretizations  based  on  our  approach  can  be  expected  to  resolve  self-sharpening 
fronts  monotonically  for  one- dimensional  reactive  transport  problems  as 
long  as  the  spatial  discretization  is  sufficiently  fine  to  allow  three  intervals 

on  a  front. 

•  The  RK2-S  tracking  approach  points  provides  a  straightforward  procedure 
for  tracking  numerical  integration  points  that  can  perform  well  for  Courant 
numbers  several  times  larger  than  one. 

•  The  RK2-FT  tracking  approach,  while  more  involved  than  the  RK2-S  al¬ 
gorithm,  allows  our  ELLAM  approaches  to  take  very  large  time  steps  and 
still  produce  accurate  solutions. 
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